{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import division, print_function\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview\n",
    "\n",
    "- http://scikit-image.org/docs/stable/api/skimage.transform.html\n",
    "- http://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.warp\n",
    "- http://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.AffineTransform (and other similar classes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image rotation from scratch\n",
    "\n",
    "The following code shows how to rotate an image using the skimage (scikit-image) library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from skimage import transform, data\n",
    "\n",
    "camera = data.camera()\n",
    "rotated = transform.rotate(camera, 30)\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ax0.imshow(camera, cmap='gray')\n",
    "ax1.imshow(rotated, cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"font-size: 150%;\">**Exercise:** Write an algorithm from scratch that will\n",
    "do the same (i.e., take an input image as an ndarray, and rotate it).</span>\n",
    "\n",
    "If you feel creative, you can also write code to magnify (zoom) the image.\n",
    "<p></p>\n",
    "You may need: http://en.wikipedia.org/wiki/Polar_coordinate_system\n",
    "<p></p>\n",
    "A (bad) solution is given below--but try it yourself before looking!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A problematic approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage import color"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rotate(image, theta):\n",
    "    theta = np.deg2rad(theta)\n",
    "    \n",
    "    height, width = image.shape[:2]\n",
    "    out = np.zeros_like(image)\n",
    "    \n",
    "    centre_x, centre_y = width / 2., height / 2.\n",
    "    \n",
    "    for x in range(width):\n",
    "        for y in range(height):\n",
    "            \n",
    "            x_c = x - centre_x\n",
    "            y_c = y - centre_y\n",
    "            \n",
    "            # Determine polar coordinate of pixel\n",
    "            radius = np.sqrt(x_c**2 + y_c**2)\n",
    "            angle = np.arctan2(y_c, x_c)\n",
    "            \n",
    "            new_angle = angle + theta\n",
    "            \n",
    "            new_x = radius * np.cos(new_angle)\n",
    "            new_y = radius * np.sin(new_angle)\n",
    "            \n",
    "            new_x = new_x + centre_x\n",
    "            new_y = new_y + centre_y\n",
    "            \n",
    "            if (new_x >= width) or (new_x < 0) or\\\n",
    "               (new_y >= height) or (new_y < 0):\n",
    "                    continue\n",
    "            else:\n",
    "                out[int(new_y), int(new_x)] = image[y, x]\n",
    "    \n",
    "    return out\n",
    "\n",
    "rotated = rotate(camera, 40)\n",
    "    \n",
    "plt.imshow(rotated, cmap='gray', interpolation='nearest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## And while we can attempt to fix the problem...\n",
    "\n",
    "...this is not an optimal approach"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Attempt at fixing the holes using a median filter\n",
    "# -- it works, sort of, but it's not the best approach.\n",
    "\n",
    "height, width = rotated.shape[:2]\n",
    "\n",
    "out = rotated.copy()\n",
    "\n",
    "for x in range(1, width - 1):\n",
    "    for y in range(1, height - 1):\n",
    "        if out[y, x] == 0:\n",
    "            out[y, x] = np.median([out[y, x-1],\n",
    "                                   out[y, x+1],\n",
    "                                   out[y+1, x],\n",
    "                                   out[y-1, x]])\n",
    "            \n",
    "plt.imshow(out, cmap='gray', interpolation='nearest');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = np.array([[4, 2], [1, 6]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(A, cmap='gray', interpolation='nearest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# For later discussion: interpolation\n",
    "\n",
    "## Bi-linear interpolation\n",
    "\n",
    "<img src=\"Bilinear_interpolation.png\" style=\"float: left;\"/>\n",
    "<div style=\"clear: both;\"/>\n",
    "\n",
    "Also see [bilinear interpolation on Wikipedia](http://en.wikipedia.org/wiki/Bilinear_interpolation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some warping experiments!\n",
    "\n",
    "## Fish-eye"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage import transform, data, io\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load face\n",
    "face = io.imread('../images/stefan.jpg')\n",
    "\n",
    "# Get the eye nicely in the middle\n",
    "face = face[:185, 15:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.imshow(face)\n",
    "plt.plot([face.shape[1]/2.], [face.shape[0]/2.], 'or', markersize=14, alpha=0.4)\n",
    "plt.axis('image');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define a transformation on the x-y coordinates\n",
    "\n",
    "def fisheye(xy):\n",
    "    center = np.mean(xy, axis=0)\n",
    "    xc, yc = (xy - center).T\n",
    "\n",
    "    # Polar coordinates\n",
    "    r = np.sqrt(xc**2 + yc**2)\n",
    "    theta = np.arctan2(yc, xc)\n",
    "\n",
    "    r = 0.8 * np.exp(r**(1/2.1) / 1.8)\n",
    "\n",
    "    return np.column_stack((r * np.cos(theta), r * np.sin(theta))) + center"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Warp and display\n",
    "\n",
    "out = transform.warp(face, fisheye)\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "ax0.imshow(face)\n",
    "ax0.set_axis_off()\n",
    "\n",
    "ax1.imshow(out)\n",
    "ax1.set_axis_off()\n",
    "\n",
    "plt.title('Knock! Knock!')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the following scripts for fun:\n",
    "\n",
    "(Open up the terminal in the \"scripts\" directory first)\n",
    "\n",
    "- **deswirl.py** (run using: ``python deswirl.py``)\n",
    "\n",
    "    In the UK, a criminal tried to hide his identity by posting\n",
    "    swirled pictures of his face online.  Here, we use the\n",
    "    Mona Lisa to illustrate what he did.  Can you restore\n",
    "    her face back to normal? (Note that you can adjust the\n",
    "    position of the red dot, as well as move the sliders.)\n",
    "    \n",
    "    \n",
    "- **clock_deblur.py**\n",
    "\n",
    "    I took a picture of a wall clock while moving the camera.  Or perhaps the clock moved.\n",
    "    Either way, now I cannot read the time!  I've implemented a deblurring\n",
    "    algorithm--can you adjust its parameters to help me pin-point\n",
    "    the time?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Here's code for a swirl transform:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage import transform\n",
    "\n",
    "def swirl(xy, center=[0, 0], strength=1, radius=100, rotation=0):\n",
    "    \"\"\"Compute the coordinate mapping for a swirl transformation.\n",
    "\n",
    "    \"\"\"\n",
    "    x, y = xy.T\n",
    "    x0, y0 = center\n",
    "    rho = np.sqrt((x - x0)**2 + (y - y0)**2)\n",
    "\n",
    "    # Ensure that the transformation decays to approximately 1/1000-th\n",
    "    # within the specified radius.\n",
    "    radius = radius / 5 * np.log(2)\n",
    "\n",
    "    theta = rotation + strength * \\\n",
    "            np.exp(-rho / radius) + \\\n",
    "            np.arctan2(y - y0, x - x0)\n",
    "\n",
    "    xy[..., 0] = x0 + rho * np.cos(theta)\n",
    "    xy[..., 1] = y0 + rho * np.sin(theta)\n",
    "\n",
    "    return xy\n",
    "\n",
    "\n",
    "h, w = face.shape[:2]\n",
    "\n",
    "parameters = {'center': [w/2., h/2.],\n",
    "              'strength': 8,\n",
    "              'radius': 90,\n",
    "              'rotation':  0}\n",
    "\n",
    "out = transform.warp(face, swirl, parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n",
    "\n",
    "ax0.imshow(face)\n",
    "ax1.imshow(out);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Can you come up with an even better distortion?\n",
    "\n",
    "## Start with this template:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def my_warp(xy):\n",
    "    x = xy[:, 0]\n",
    "    y = xy[:, 1]\n",
    "    \n",
    "    x = x + 1.5 * np.sin(y / 3)\n",
    "    \n",
    "    return np.hstack((x, y))\n",
    "\n",
    "image = plt.imread('../images/stefan.jpg')\n",
    "out = transform.warp(image, my_warp)\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax0.imshow(image)\n",
    "ax1.imshow(out);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Composing Transformations\n",
    "\n",
    "scikit-image allows you to compose several transformations.  For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage import data\n",
    "\n",
    "cat = data.chelsea()\n",
    "horizontal_shift = transform.SimilarityTransform(translation=[20, 0])\n",
    "\n",
    "multiple_shifts = horizontal_shift + horizontal_shift + horizontal_shift\n",
    "\n",
    "f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "ax0.imshow(cat)\n",
    "ax1.imshow(transform.warp(cat, horizontal_shift.inverse))    # Note the inverse!\n",
    "ax2.imshow(transform.warp(cat, multiple_shifts.inverse));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `transform` module allows us to rotate images.  The inner workings is something like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def my_rotate(image, angle):\n",
    "    rotation_tf = transform.SimilarityTransform(rotation=np.deg2rad(angle))\n",
    "    return transform.warp(image, rotation_tf.inverse)\n",
    "\n",
    "plt.imshow(my_rotate(cat, 30))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this rotates the cat around the origin (top-left).\n",
    "\n",
    "**Can you modify `my_rotate` to rotate the image around the center?**\n",
    "\n",
    "*Hint:*\n",
    "\n",
    "1. Shift the image (see above) so that the center of the image lies at (0, 0)\n",
    "2. Rotate the image\n",
    "3. Shift the image back---the opposite of what you did in step 1\n",
    "\n",
    "All of this can be achieved by composing transformations and calling `warp` once."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced challenge: rectifying an image\n",
    "\n",
    "<img src=\"../images/chapel_floor.png\" style=\"float: left;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know the above tiles are laid out in a square--can you transform\n",
    "the image so that the tiles are displayed as if you were viewing them from above?\n",
    "\n",
    "The centre-points of the corner circles are, given as (row, column) coordinates:\n",
    "\n",
    "```\n",
    "(72, 129) -- top left\n",
    "(76, 302) -- top right\n",
    "(185, 90) -- bottom left\n",
    "(193, 326) -- bottom right\n",
    "```\n",
    "\n",
    "Hint: there is a linear transformation matrix, $H$, such that\n",
    "\n",
    "$H \\mathbf{x} = \\mathbf{x}'$\n",
    "\n",
    "where $\\mathbf{x}$ is the *homogeneous* coordinate in the original image and\n",
    "$\\mathbf{x}'$ is the *homogeneous* coordinate in the rectified image (with *homogeneous*\n",
    "we simply mean that we add an extra 1 at the end, e.g. (72, 129) becomes (72, 129, 1).\n",
    "The values for $\\mathbf{x}$ and their new values, $\\mathbf{x}'$,\n",
    "are therefore:\n",
    "\n",
    "```\n",
    "x = (72, 129, 1), x' = (0, 0, 1)\n",
    "x = (76, 302, 1), x' = (0, 400, 1)\n",
    "x = (185, 90, 1), x' = (400, 0, 1)\n",
    "x = (193, 326, 1) x' = (400, 400, 1)\n",
    "```\n",
    "\n",
    "(You can choose any output size you like--I chose $400 \\times 400$)\n",
    "\n",
    "Why do we need homogeneous coordinates?  It allows us to have *translation* as part of H:\n",
    "\n",
    "$\n",
    "\\left[\\begin{array}{ccc}\n",
    "H_{00} & H_{01} & H_{02}\\\\\n",
    "H_{10} & H_{11} & H_{12}\\\\\n",
    "H_{20} & H_{21} & 1\n",
    "\\end{array}\\right]\\left[\\begin{array}{c}\n",
    "x\\\\\n",
    "y\\\\\n",
    "1\n",
    "\\end{array}\\right]=\\left[\\begin{array}{c}\n",
    "H_{00}x+H_{01}y+H_{02}\\\\\n",
    "H_{10}x+H_{11}y+H_{12}\\\\\n",
    "H_{20}x+H_{21}y+H_{22}\n",
    "\\end{array}\\right]\n",
    "$\n",
    "\n",
    "Note that each element of the output coordinate is of the form $ax + by + c$.  Without the 1 in the last position of the coordinate, there would have been no $+ c$ and therefore no translation!\n",
    "\n",
    "The question on how to determine $H$ is left for another day.  If you are curious, \n",
    "the [answer can be found here](homography.pdf).\n",
    "\n",
    "In the meantime, I provide some code to calculate $H$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from skimage.transform import estimate_transform\n",
    "\n",
    "source = np.array([(129, 72),\n",
    "                   (302, 76),\n",
    "                   (90, 185),\n",
    "                   (326, 193)])\n",
    "\n",
    "target = np.array([[0, 0],\n",
    "                   [400, 0],\n",
    "                   [0, 400],\n",
    "                   [400, 400]])\n",
    "\n",
    "tf = estimate_transform('projective', source, target)\n",
    "H = tf.params\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the code in the cell above, you can compute the target coordinate of any position in the original image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Verify that the top left corner maps to (0, 0)\n",
    "\n",
    "x = np.array([[129, 72, 1]])\n",
    "\n",
    "z = np.dot(H, x.T)\n",
    "z /= z[2]\n",
    "\n",
    "print(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here's a template solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def rectify(xy):\n",
    "    x = xy[:, 0]\n",
    "    y = xy[:, 1]\n",
    "    \n",
    "    # We need to provide the backward mapping, from the target\n",
    "    # image to the source image.\n",
    "    HH = np.linalg.inv(H)\n",
    "    \n",
    "    # You must fill in your code here to take\n",
    "    # the matrix HH (given above) and to transform\n",
    "    # each coordinate to its new position.\n",
    "    # \n",
    "    # Hint: handy functions are\n",
    "    #\n",
    "    # - np.dot (matrix multiplication)\n",
    "    # - np.ones_like (make an array of ones the same shape as another array)\n",
    "    # - np.column_stack\n",
    "    # - A.T -- type .T after a matrix to transpose it\n",
    "    # - x.reshape -- reshapes the array x\n",
    "    \n",
    "    # ... your code\n",
    "    # ... your code\n",
    "    # ... your code\n",
    "    # ... your code\n",
    "    \n",
    "    return ...\n",
    "\n",
    "    \n",
    "image = plt.imread('images/chapel_floor.png')\n",
    "out = transform.warp(image, rectify, output_shape=(400, 400))\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax0.imshow(image)\n",
    "ax1.imshow(out)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"height: 100px;\">&nbsp;</div>\n",
    "\n",
    "<div>\n",
    "The solution to the above problem is provided as [solutions/tile_rectify.py](solutions/tile_rectify.py).  Only look at it after you've attempted the problem yourself!\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# For more fun examples see http://scikit-image.org/docs/dev/auto_examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1+"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
